o 

(N 
O 

Q 



U 

o 

o3 



Scenario Trees and Policy Selection for Multistage 
Stochastic Programming using Machine Learning 

Boris Defourny 

Department of Operations Research and Financial Engineering, Princeton University, 
Princeton, New Jersey 08544, USA, defourny@princeton.edu 

Damien Ernst, Louis Wehenkel 

Department of Electrical Engineering and Computer Science, University of Liege, 
Grande Traverse 10, 4000 Liege, Belgium, {dernst@ulg.ac.be, L.Wehenkel@ulg.ac.be} 

We propose a hybrid algorithmic strategy for complex stochastic optimization problems, 
which combines the use of scenario trees from multistage stochastic programming with ma- 
chine learning techniques for learning a policy in the form of a statistical model, in the context 
of constrained vector-valued decisions. Such a policy allows one to run out-of-sample simu- 
lations over a large number of independent scenarios, and obtain a signal on the quality of 
the approximation scheme used to solve the multistage stochastic program. We propose to 
apply this fast simulation technique to choose the best tree from a set of scenario trees. A 
solution scheme is introduced, where several scenario trees with random branching structure 
are solved in parallel, and where the tree from which the best policy for the true problem 
could be learned is ultimately retained. Numerical tests show that excellent trade-offs can 
be achieved between run times and solution quality. 
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1 Introduction 

Stochastic optimization using scenario trees has proven to be a powerful algorithmic strategy, 
but ha s suffered from the rapid growth in the size of sc enario trees as the number of stages 
grows (IBirge and Louveauxl . 119971 : IShapiro et all l2009l ). A number of authors have under- 
taken research to limit the size of the scenario tree, but problem size s t ill grows exponentially 
with the number of stages (IFrauendorfe r. 1996 ; iDupacova et al.l . |2000| ; iHpyland and Wallace! . 
200ll ; iPennanenl . 120091 ; iHeitsch and Romischl . 120091 ) . As a result, most authors either s everely 
limit the number of decision stages or sharply lim i t the number of scenarios p er stage (IBirge , 
19971 ; IWallace and Ziembal . 120051 ; iDempster et al.l . 120081 ; iKallrath et al.l . 120091 ) . Such approx- 
imations make it possible to optimize first-stage decisions with a stochastic look-ahead, but 
without tight guarantees on the value of the computed decisions for the true multistage 
problem (as a matter of fact, bounding techniques also tend to break down on problems 
with many stages). 

Some authors hav e proposed to asses s the quality of sc e nario - tree based methods by out - 
of-sample validation (jKouwenbergl . l200ll ; IChiralaksanakull . l2003t iHilli and Pennanenl . |2008|). 
The validation scheme consists of solving the multistage program posed on a scenario tree 
spanning the planning horizon T, implementing the decision relative to time step 1, sampling 
the realization of the stochastic process at time 1, updating the conditional distributions of 
the stochastic process from time 2 to time T, rebuilding a scenario tree spanning time peri- 
ods 2 to T, solving the new multistage program over the remaining horizon (with previously 
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implemented decisions fixed to their value), and continuing this process until the last deci- 
sion at time T has been found. The resulting sequence of decisions is then valued according 
to the true objective function. Averaging the result of this procedure repeated over many 
independent scenarios drawn according to the true distributions of the problem produces an 
unbiased estimate of the expected value of the solution for the true problem. Unfortunately, 
such simulations are very demanding computationally. Moreover, the variance of the em- 
pirical estimate is likely to be larger for problems with many stages, calling for even more 
simulations in that case. As a result, running times restrict the use of this technique to 
relatively simple optimization problems and simple scenario-tree updating schemes. 

In this paper, we propose a hybrid approach that combines scenario trees with the es- 
timation of statistical models for returning a decision given a state. One could call these 
statistical models policy function approximations (PFA). We solve an optimization problem 
on a scenario tree to obtain optimal decisions (given the tree), and then we use the decisions 
at all stages in the t ree to fit the policy function approximations using methods based on 
supervised learning (IHastie et all 120091 ). We repeat this exercise for different samples of 
scenario trees, producing a family of policy function approximations. Each policy function 
approximation is then tested on a fresh set of samples to determine the best policy function 
approximation. 

Machine learning methods have often been a pplied to stochastic optimization, primarily 

in the context of approximating a value func t ion (IBerstsekas and Tsitsiklisl. 1 1 996c ISutton and Bartol . 



19981 : iBertsekasl . l2005t iBusoniu et all 1201(1 ISzepesvaril . 12010 : iPowelll . 1201 if ). The statistical 



estimation of policies has also been widely stud ied in the reinforcement learning community , 



often using the term "actor-critic" methods (ISutton and Bartol . Il998t iPeters and Schaal 



20081 ). Such methods are popular in computer science for discrete action spaces, and in con- 



trol theory for low- dimensional but unconstrained control problems. Our method is designed 
for higher- dimensional constrained optimization problems. Since we cannot capture compli- 
cated constraints in the policy function approximations, we solve a constrained optimization 
problem for finding the best feasible solution that minimizes the deviation from the PFA. 

We note that some authors have also proposed to derive a policy from a scenario tree 
by applying to a new s c enario the decision optimize d for the closest scenario in the tree 



(IThenie and Viall . l2008t iKiichler and Vigerskd . |2010J); their strateg y could be viewed as 



a form of apprenticeship learn i ng by nearest neighbor regression (lAbbeel and Ngi . 12004 



Syed et all . 120081 ; ICoates et all . 120081 ). However, the use of machine learning along with 



a valid model selection procedure is quite new in the context of stochastic programming, 
while the need for methods able to discover automatically good decision rules had been 
recognized as an important research direction for addressing complex multistage stochas- 
tic programmin g problems (IMulvey and Kiml . 1201 ll ) and for bounding approximation errors 
flShapirol - lioblh . 

The machine learning approach makes it possible to quickly perform out-of-sample eval- 
uations of the policy function approximations created using each scenario tree. The result 
is that it is now practical to compute safe statistical guarantees when using a scenario tree 
approximation scheme. Building on this ability, we propose to revisit the fundamental prob- 
lem of generating the scenario tree from which a first-stage decision to be implemented is 
optimized. We are particularly interested in working with trees that have a sparse branching 
structure for representing uncertainty over a long planning horizon using a limited number 
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of scenarios. Optimization over long planning horizons is especially relevant for exploiting 
a resource in limited quantity such as water in a reservoir, or an electricity swing option, 
when the price of the resource is stochastic. Small trees are also a pragmatic choice when 
the number of scenarios is very limited by the complexity and dimension of the problem, for 
instance in stochastic unit commitment problems for electricity generation scheduling. 

In this paper, we consider a randomized algorithm for generating small scenario trees 
over long horizons, that uses a branching process for generating the branching structure. 
We illustrate the solution approach on a set of problems over a long planning horizon. A 
small fraction of the PFAs learned from the random scenario trees turn out to perform very 
well in out-of-sample simulations. We need not ask more from the scenario-tree generation 
algorithm, in sharp contrast with solutions approaches based on the optimization of a single 
scenario tree. 

Our approach to scenario tr ee generation can be seen as an extension of the Stochastic 
Approximation Method (SAA) (IShapiro et all 120091 ). where in addition to drawing scenarios 
randomly, the branching structure of the tree is also drawn randomly. However, as the 
detection of good structures is left to the out-of-sample validation, the way of thinking about 
the scenario trees is radically changed: in our solution approach, the best-case behavior of 
the scenario tree generation algorithm can be exploited. Our approach could be contrasted 
to other solution schemes based on multiple tr ees (possibly each reduced to one scenario ) 
used inside averaging or aggregation procedures (IMak et al.l . Il999t iNesterov and Viall . 120081 ) , 
perhaps most notably in iRockafellar and Wetsl ( 119911 ) which inspired the title of this paper. 

Our paper makes the following contributions. We introduce the hybrid policy structure 
along with a model selection procedure for selecting a best policy function approximation, 
given a scenario tree solution and the specification of the probabilit y distributions and con - 
straints of the true problem. This idea was originally presented in iDefourny et al.l ( 120091 ) , 
where complexity estimates were given (but without working algorithms). We identify sta- 
tistical models amenable to fast simulation in the context of convex multistage stochastic 
programming problems, so as to be able to quickly generate and simulate feasible decisions on 
a large test sample. We conduct numerical tests on a multistage stochastic problem to assess 
the sensitivity of the approach to various choices, and its benefit in terms of running time 
for obtaining a solution with a statistical performance guarantee. We propose a novel way to 
view the problem of constructing the scenario trees, which has the potential to scale better 
with the number of decision stages without imposing stringent conditions on the problem 
structure. We report successful numerical experiments obtained with this technique. 

The remainder of this paper is organized as follows. Section [2] explains the solution 
approach. Section [3] formalizes the description and gives algorithmic procedures for its 
implementation. Section H] investigates the method numerically on a test problem. Section [5] 
describes the proposed application of machine learning techniques to scenario tree selection, 
reports numerical results, and Section [6] concludes. 



2 Principle 

Let us explain the principle of our solution approach on a stylized example. Consider decision 
stages numbered from t = 1 to t = T with T large, say T = 50. Assume that the decision at 
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stage t is a scalar x t G X = [0, 1] that can be adapted to the history h t = (£1, £ 2 , . . . , £ t ) = 
£m of a random process £ t . Assume that £1 is constant and (for concreteness) that £ 4 for 
t > 1 follows a standard normal distribution. We can write the multistage problem as 
V : minE/(xi, £ 2 , x 2 (/i 2 ), £t, ^t(^t)), assuming that one observes £ t and then takes 
decision x t . Let us denote by v the optimal value of the problem, assumed to be finite. 

Assume first that a scenario tree T is given to us. Each node j of the tree represents 
an information state h(j) = (£i, . . . ,£*(>•)) where is the stage determined by the depth 
of node j. Transition probabilities are associated to the arcs between successor nodes. Each 
leaf node k of the tree determines, by its path from the root node, a particular scenario 
£ fc = £f, • • • j £f0- The probability of the scenario, written p k , is obtained by multiplying 
the transition probabilities of the arcs of the path. Thus, the probability of reaching an 
information state h(j) is the sum of the probabilities p h of the scenarios passing through 
that node. Note that most nodes can have only one successor node, since otherwise the 
number of scenarios would be an astronomical number on this problem with so many stages. 

On the scenario tree, we can formulate a math program where optimization variables 
are associated to the nodes of the tree. For each scenario k, we associate optimization 
variables x\, . . . , x\ to the nodes on the path from the root to the leaf k, and we enforce 
identity among the optimization variables that share a common node in the tree (non- 
anticipativity constraints). We solve V(T) : min^2 k p k f(x , l, £f,2^> . . . , £f., xff) subject to 
the non-anticipativity constraints. Let v(T) and x*(T) denote its optimal value and optimal 
first-stage decision. 

At this stage, the regret of implementing x^(T), that is, 

minE/(;ri, ^2,^2(^2), £tj Xr(hr)) — v subject to x\ = x\{T) , 

is unknown: the tree represents £ t given . . . , £t_i) by a single realization at most of the 
nodes, so the regret depends on the value of the stochastic solution for the subproblem at 
each node. 

View now Xt as a mapping from the history h t G M* to a decision in X = [0,1]. The 
optimal solution on the scenario tree provides examples of input-output pairs (h^,x k ). By 
machine learning, we can infer (learn) a statistical model for each mapping x t . In the 
paper this idea will be generalized to vector-valued decisions and convex feasibility sets by 
building one model per coordinate and restoring the structure of the vector by solving a small 
math program. For simplicity, we illustrate the idea in the scalar case and with parametric 
regression. By a set of points in the input-output space R* x [0, 1], we can fit a constant- 
valued function, a linear function, a quadratic function, and so forth (Figured]). We can 
build statistical models for xi,. . . ,Xt and view them as a single model it = (71*1, . . . , vr^) for a 
policy from stage 1 to T, where 7Ti = x± is equal to x\(T). Among a set of statistical models 
tt u , where v indexes the model, we cannot know in advance the model that works best for 
the true problem. However, we can sample N new independent scenarios (£1, • • • > £t) 
according to the true distributions (in this example, standard normal distributions) and 
guarantee, if N is sufficiently large, that 

v < min Eli /K Ta(£i, £2), •••>£?, ^(£1, £ 2 n , • • • , £?)) + 0{N~^) } , 
where the right-hand side is simple to evaluate due to the nature of the models ix v t . 



4 



dataset [0, 1] 




Figure 1: Possible statistical models for a recourse function x t (h t ) where h t = . . . , £ t ), 
from data extracted from a scenario tree: Illustration with parametric regression. The 
superscript indicates the scenario. 

The next step is now to have the whole process repeated for a large number of different 
trees, where each tree has potentially a different first-stage decision. That is, if we have an 
algorithm A that maps a vector of parameters u m to a scenario tree T m , we can generate 
M trees by selecting or sampling u m , and from each tree T m learn statistical models ir u ' m 
from the solution to V(T m ). Using a large sample of N scenarios, we can guarantee, if N is 
sufficiently large, that 

v < min Eli ®, <' m (6, 8), • • • , ©, 7# m (£i, . . . , ©)) + 0(N-W)\ . 

We can then implement the decision 7r™ of the tree that attains the minimum, or even use 
the best model 7i m ' u over the whole horizon. Note that an unbiased estimate of the expected 
value of the selected model ii m ' u on the true problem is obtained by simulating the model 
again on a new independent test sample of scenarios. 

For trees over a large planning horizon, we propose to consider algorithms A capable of 
generating different branching structures given a target number of scenarios, so as to keep 
the complexity of V(l~ m ) under control. 

In the remainder of the paper, we formalize this approach, using nonparametric statis- 
tical models, and develop methods for dealing with vector-valued decisions and feasibility 
constraints. A first set of numerical experiments compares the bounds obtained with the 
machine learning approach to the ideal bound computed by solving a multistage program 
at each stage for each out-of-sample scenario (on a problem over 4 stages). A second set of 
numerical experiments evaluates the approach on a problem over 50 stages, in combination 
with a simple implementation of A based on a branching process that modifies itself its 
branching probabilities to control the total expected number of scenarios. 

As a last remark, we note that the optimal values v{T m ) of the programs V{T m ) are a 
poor indicator of the relative true quality of the corresponding solutions 7r™. Our solution 
approach does not use v(T m ), and in fact would still be valid with trees constructed by 
tweaking the true distributions. This could be useful for instance to inflate the probability 
of generating extreme event scenarios. 
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3 Mathematical Formalization 



This section formalizes the supervised learning approach proposed in this paper. After 
summarizing the notations in Section l37Lj Section l3~2l states the generic form of the multistage 
stochastic program under consideration. Section I3U1 gives the generic form of a scenario-tree 
approximation. Section 13.41 describes the datasets extracted from an optimal solution to the 
scenario-tree approximation. Section 13.51 describes a Gaussian Process regression method 
that infers from the datasets a statistical model of the mapping from information states to 
recourse decisions. Section 13.61 describes the optimization-based method that exploits the 
statistical model to output a feasible decision at each stage, given the current information 
state. Section 13.71 describes the simulation-based procedure for selecting a best statistical 
model in combination with the feasibility restoration procedure. 



3.1 Notations 



We follow notations proposed in IShapiro et al.l (120091 ). Equation 3.3, where the random 
process starts with a random variable £ x that has a single trivial value. 

t : stage index, running from 1 to T. 

x t G M. nt : decision at stage t. 

£ t : random vector observed just before taking decision xt- 

6[t] = • • ,6): the random vectors up to stage t, where £i has a single trivial value. 
That is, the history of the random process observed at stage t, which represents the 
information state at stage t. 

/t(-): cost function at stage t, depending on x t and £[ f ]. 

x u x 2{-)i ■ ■ ■ , x t{-)' decision policy from stage 1 to T, where the value of x t is uniquely 
determined by the realization of £[ t j. For defmiteness, the domain of x t (-) must comprise 
the support of £w. 

X t : feasibility set for x t , in the sense that x t must be valued in X t . The set X t may be 
expressed using a convenient representation of the information state at time t. 

(a, b) : inner product a T b between two vectors a, b of same dimension. 



3.2 True Problem 

The multistage stochastic program under consideration is called the true problem. The true 
problem is written in abstract form as 

min^,...^ E + /aMfp]), 6) +■■■ + /r(Er(f[r|),6r)] ^ 

subject to Xi G X u z t (£[ t ]) G rY t (x t _i(C[t-i]), t = 2,...,T. 
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We also keep in mind Theorem 2.4 in iHeitsch and Romischl (120111 ). which establishes 



assumptions for ensuring that perturbations of the e-optimal solution sets of a program such 
as ([I]) remains bounded under small perturbations of the distributions for £[ t j. A simple 
particular case of the theorem works with the following assumptions: 

1. The multistage stochastic program is linear, with 

fi(xi) = (Ci,Xi), 

/t (a* (£[*]) ,6) = (<H(tt),x t (€[t])), 

X x = {xi G X\ : A^qXx = rjt}, 

Xt(x t -i[Z[ t -i]),€t) = {%t € X t : A tfi x t + A t>1 x t -i = ?7i(£t)} 

for some fixed vectors ci, rjx, and for some vectors c t , r] t depending affinely on £ f , 
t = 2, . . . , T, and for some nonempty, bounded, fixed sets X t and some fixed matrices 
A tfi ,t = l,...,T. 

2. £i and #t (£[*]) have finite second moments. 



3. IHeitsch and Romischl (120111 ) also specify a perturbation-robust version of the relatively 
complete recourse assumption. Under small bounded perturbations of the conditional 
distributions for £ t given £[t-i]> it must still be the case that the relatively complete 
recourse assumption holds. 



4. IHeitsch and Romischl (120111 ) also give a perturbation-robust version of assumptions 



ensuring that the optimal value of the perturbed problem is attained. 

3.3 Approximation of the True Problem 

The scenario-tree approximation for (TTT) is written in abstract form as 



subject to x\ e X x , x\ e X^x^,^), t = 2, . . . , T, k = 1, . . . , K, (2) 
x\ = x[ for all k, I such that = £L, t — 1, . . . , T. 

Here x\, . . . , x\ denote decision vectors relative to scenario k for k = 1, . . . , K, £fc is the 
history of random vectors up to stage t for scenario k, p k is the probability of scenario k 
with Y^k=iP k = 1> an< i t ne additional constraints are the nonanticipativity constraints. 

The discrepancy between the optimal value, solution set of (j2J) and the optimal value, 
solution set of ([1]) depends on the number K of scenarios, the branching structure of the 
tree, and the values p k , £|L. Conditions for the epi- convergence of to (pQ) (con vergence of 



the op timal value and the optimal set for the first-stage decision) are studied in iPennanen 



(120051 ). building on the work of many others. 
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3.4 Extraction of State-Decision Pairs 



Let us define datasets of state-decision pairs extracted from an optimal solution to (j2J): 

A = *f)}£=i> t = l,...,T. (3) 

Here x±*, . . . , xij? denote the optimized decisions relative to scenario k. Note that the dataset 
D\ always contains the same pair (£i,x*) representing the constant first-stage decision. The 
dataset Dt contains distinct pairs. The datasets D t for t — 2, . . . , T — 1 have some duplicate 
entries resulting from the branching structure of the scenario tree, assuming all necessary 
non-anticipativity constraints have been formulated and §Z§ has been solved to optimality. 
We let \D t \ < K denote the number of distinct pairs in D t . 

Clearly, the pairs in the datasets are not independent, identically distributed (i.i.d.) 
samples from a fixed but unknown distribution, since the decisions x\ are optimized jointly. 
Recall also that the optimized decisions are not the optimal decisions for the true problem. 

Depending on the machine learning strategy, it can be beneficial to represent the in- 
formation state and decision by minimal representations, or expand the representation to 
express additional features. 



3.5 Inference of the Statistical Model for Making Decisions 



Several statistical models are possible, but in the present context we find it particularly 
attractive to consider nonparametric models. We work with Gaussian processes, which 
are relatively independent of the way datasets have been generated (in particular the i.i.d. 
assumption). With Gaussian processes, prior distributions ove r spaces of functions can 
be d efined and then updated given input-output observations (IRasmussen and Williams! . 
20061 ). The prior over decisions is determined by the choice of covariance functions Cf (•, •) 
parameterized by a vector 9 and by the choice of mean functions mt(-), for t — 2, . . . , T. For 
a decision at stage t described by n t scalar variables, the function m t { ) is vector- valued with 
n t coordinates written m t i (-), i = 1, . . . , n t . 

We define the following short-hand notations: 



<i(A) 

m t A D t) 

C?(A,A) 



X 



X 



\D t \* 



t.i 



m tA[ t ]) ■■ 



Dt|xl 

) 

T 



m (£ {Dtl 



r<6(£l t\D 

°i ^[t]'S[i] 



Dt\ C l 



\D t 



\D t \xl 



\D t \xl 



\Dt\x\D t \ 



(4) 
(5) 
(6) 

(7) 



Examples of ker nels, analysis of their prop e rties, and vali d ways to combine several k ernels , 
can be found in Rasmussen and Williams fl2006h and in ISteinwart and Christmanl ( 120081 ). 
The most important property is that Cf(D t ,D t ) must be positive semidefinite for any D t . 
Intuitively, the kernel induces a similarity measure in the input space. 
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We will also assume that the samples x k * in the datasets ([3]) are noisy observations 
following the (nonparametric) statistical model 

x k t * = ^\ t] ) + w k t (8) 

where n t denotes an unknown but fixed optimal policy at stage t for the true problem, is 
generated according to the distributions specified in the true problem, and w\ is a zero-mean 
Gaussian noise of variance cx^. 

The statistical model jSJ) is an approximation that can be motivated by an asymptotic 
argument. With a scenario-tree approximation method such that d2J) epi-converges to (1T1) 



(jPennanenl . 120051 ) . the decisions x k * optimal for d2J) tend to decisions nt(& t ]) optimal for §Q 



as the number of scenarios K and the number of branchings in the tree tend to infinity. 
Then, with an infinite number of examples covering the state-decision space, we can select 
a radial basis kernel with a bandwidth tending to 0, and thus construct a nonparametric 
model for ir t that replicates the examples arbitrarily well. 

Given we treat all the coordinates of a stage- 1 decision independently. Under the 
updated distribution over policies, the predictive distribution for the unknown tt^ evaluated 
at £m is Gaussian with mean and variance given respectively by 

Aye M )=m t)i (£ [t] ) + C?(^^ , (9) 

= c?(£ w ,£w) - c?(A,£w) T [c?(A, A) + (T^cfCA.ew) , (io) 

with / denoting the identity matrix in Ml D *l x l Dt l. 

If we want to replicate closely the decisions optimal for fl2]), a small noise variance <r^ 
should be chosen, that will simply act as a numerical regularizer in the inversion of the 
matrix Cf(D t ,D t ). If we have some confidence in the prior, or equivalently some mistrust 
in the approximation of (pQ) by a larger noise variance could be chosen. A larger noise 
variance reduces the weight of the updates made to the decisions m tt i(£[t]) determined by the 
prior. 

For the mean functions, a common choice in the machine learning literature is to set 
m t = 0. This reflects a noninformative prior. Another option would be to solve first a 
deterministic model (typically the multistage problem on a single nominal scenario), extract 
the optimized decisions x t , and set m t = x t . 

The choice of the covariance functions Cf (along with their parameters 9) affects the 
"regularity" of the decision policies that are generated from the updated distribution. There 
may exist policies optimal for the true problem that exhibit kinks and discontinuities as 
functions of £[ t ], but usually (that is, with mean functions set to 0) it is not possible to 
obtain discontinuous functions from Gaussian processes. In the present paper, the choice of 
the covariance function and of their parameters 9 is incorporated to a general model selection 
procedure, based on the simulation of the decision policy on the true problem (see §3.7p . 

3.6 Inference of a Feasible Decision Policy 

In a Bayesian framework, the estimate of the decision a^(£[t]) G M. nt is not described by a 
single vector, but by a predictive distribution. 
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Figure 2: Restoring the feasibility of a prediction \ t for x t G X t by solving (fTTj) . where a* 
denotes the corresponding optimal value. 

We define in this section simple selection procedures that will output a single feasible 
decision x t (£[t])- 

We know that under the Gaussian process model, the density of the predictive distribution 
is the density of a multivariate Gaussian with mean \f(£,[t]) = [^ti(C[t])j ■ ■ ■ ■> nt(^[t})V an d 
covariance matrix A? (£[*]) = diagfAj^ffl), . . . , Aj nt (^ t] )), using Q, (JTOj). 

We can select a single feasible decision ) t° be implemented by maximizing the 
log-likelihood of the density, subject to feasibility constraints: 

z t (£ M )= argmin, t (z t - A?(£ [t] )) T [A? (^i)]"^* - A?(£ [t] )) (n) 
subject to x t E X t (x t -i(Cit-i]),^t) ■ 

The program (iTTj) is essentially the implementation of a projection operator on the feasibility 
set ^(xt-i^ft-i]), £t), applied to the conditional mean A? (£[$]): see FigureEl 

Solving (fTTj) after the evaluation of its parameters induces a feasible decision policy. 

Another, faster option is first to select the mean Af(£w) and then to correct it with 
some fast heuristic for restoring its feasibility. The heuristic could have a small number 
of parameters 9. The heuristic is interpreted as a part of the decision maker's prior on 
near-optimal decision policies. 



3.7 Model Selection 

Ultimately, a decision policy should be selected not for its ability to explain the decisions of 
the scenario-tree approximation, but for its ability to output decisions leading to the best 
possible expected performance on the true problem. 

The performance of any feasible policy x = (xi, . . . , xt) is actually the value of x for the 
true multistage problem, written v(x). It can be estimated by Monte Carlo simulation over 
a large test sample of i.i.d. scenarios £#i = . . . 1 < j < N, independent of the 

scenarios of the scenario tree. The estimator of v(x) on the test sample TS = j^*- 7 }^ is a 
sample average approximation (SAA) estimator: 



v TS (x) y 



1 N 



'12) 



10 



If v(x) is finite, the estimate (11 2p is an unbiased estimator of the value of x on the true prob- 
lem, by the strong law of large numbers. Moreover, if the objective of the true problem under 
the policy x has a finite second moment, by the central limit theorem (fl2j) is approximately 
normally distributed, with a variance that can be estimated by 

° 2 (v TS (x)) = £^[/i(*(**f)) + • • • + hmf^Up) ~ v TS (x)} 2 ) . 

Hence, one can guarantee with an approximate level of confidence (1 — a) that 

v(x) < v TS (x) + z a/2 a(v TS (i)) (13) 

with z a /2 = $ -1 (l — a/2) and denoting the invers e cumulative distribution function of 



the standard normal distribution (IShapiro et all 120091 ). 



The right-hand side of ({TBI) is a statistical performance guarantee on the true problem. 
Ranking different policies derived from different priors is possible on the basis of (I13j) . al- 
though more efficient ranking and selection techniques could be employed to eliminate more 
rapidly bad policies (we are thankful to Alexander Shapiro for this suggestion). 



4 Numerical Test 

We investigate the proposed methodology numerically according to three main factors of 
variation in its implementation: (i) the size of the scenario tree used for approximating the 
true problem, relatively to the size that should be used to solve the multistage problem accu- 
rately; (ii) the choice of the covariance function of the Gaussian processes, that determines 
how the decisions extracted from the scenario tree are generalized to new information states; 
(iii) the choice of the feasibility restoration procedure, which plays a role if the predicted 
decisions are not feasible. 

The decision policies derived from these experiments are evaluated according to two 
criteria: (i) the quality of the decision policy, relatively to the best performance attainable 
for the true problem; (ii) the computational complexity of simulating the policy. 

The experiments are imp lemented in Matlab and the programs are solved with cvx 
(iGrant and Bovdl . 120081 . l2009h . 



11 



4.1 Test Problem 



In the spirit of a stylized application presented in IShapiro et al.l ( 120091 ). Section 1.3.3, we 
consider a four-stage assembly product problem: 



min XuX2:X3:X4 E [(ci,Xx) + J2 t=2 (c t , %t(€[t] y 
subject to x x G X t , a?t(f[t]) e -**(zt-i(f[t-i])> * = 2 ;3,4, 

A!i = {xi G M 12 : ar^ > 0}, 

X 2 (xi) = {x 2 = (q 2 ,Y 2 ) G M 8 x M 12xS : 

< Y 2ji j, Y,j A 2,ij < q2,i,Y 2tij > 0} (14) 

* 3 M£ [2 ])) = {x 3 = (q 3 ,Y 3 ) G R 5 x M 8x5 : 

X i (x 3 (£ [3] ),^ ] ) = {x i em. 5 : 

Xi,i < Vi(€[4]) ■= max{0, (6i,C[4])}, < x 4 ,< < ?3,i(C[3])} 

with £1 = 1, and ^2,^3,^4 i-i-d. random variables each following a standard normal distribu- 
tion. The problem data is given in appendix (Appendix 1X1). 

In this resource allocation problem, £ represents observable factors that contribute to 
demands ?7i(£[4]) fully revealed at the last stage. The decisions q tj i represent output quantities 
of a component of type (t,i). Decisions Y tji j represent the quantity of component (t — 
allocated to the production of type (t,j), in a proportion fixed by the composition coefficient 
At,ij > 0. The cost coefficients c t have nonnegative components at stages 1, 2, 3, but 
are negative at stage 4, for representing revenue drawn from selling the 5 end-products in 
quantity at most equal to the demand or to the inventories. 

The size of the test problem has been fit to the numerical experiments to be conducted. 
For benchmarking, we simulate pure multistage stochastic programming decision processes, 
that work by instantiating and solving a new version of ( fl4|) at each decision stage, over the 
remaining horizon, with the previous decisions fixed at their implemented value. Obtaining 
the benchmark values takes many hours of computation on a single processor (Section I4.4p . 
but the simulations could run easily in parallel once the common first-stage decision has 
been computed. 

The numerical parameters of the test problem (Appendix have been chosen by select- 
ing, among randomly generated sets of parameters, a set of paramete rs that "maximizes 1 



the v alue of the multistage model. The value of the multistage model (jHuang and Ahmed 



20091 ) is the difference (in absolute value) between the optimal value the multistage model 
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(|T4|) and the optimal value of the corresponding two-stage model, namely 



min^^g^ Y?t=i{°ti x t) + E [( c 4, ^4(C[4]))] 

subject to (x 1 ,x 2 ,x 3 ) G X 1 x A^Zi) x ^(a^), x 4 (^[ 4 ]) G ^4(^3, £[4]), 

= {xi G K 12 : x M > 0}, 
X 2 { Xl ) = { x 2 = (q 2 , Y 2 ) Gl 8 x M 12x8 : 

A 2 ,ijq 2 j < Y 2ji j, V ; .1l' ,j < g 2ji , y 2 ,ij > 0} (15) 

Xs(x 2 ) = { x 3 = (q 3 ,Y 3 )eR 5 xR s » 5 : 

;t 4 (x 3 , £[4]) = { X 4 GM 5 : 

S4,i < ?7i(£[4]) := max{0, (&<,f[4])}, < Z4,<(£[4]) < 9s,*}- 

The two-stage model f lT5|) does not exploit the opportunity of adapting the production 
plan to intermediate observations available before the demand is fully revealed. 

4.2 Studied Variants in the Learned Policies 

The policy function approximation to be learned is made of two components, the statistical 
model (Gaussian processes) and the feasibility restoration procedure. This section describes 
the variants that we have considered for the tests. 

4.2.1 Covariance Functions 

We report results obtained with covariance functions of the form 

C?(£fe,$|) = exp { - [Et =1 {9(e) ~ 9(C)) 2 ] I (2 Q 2 ) } (16) 

for two choices of the function g(-) : R — > R, which must be continuous and one-to-one. 
In the first variant, g(-) is reduced to the identity function. Hence ffTol) is a radial basis 
function and 9 > is the bandwidth parameter. In the second variant, g(-) = $(-), the 
cumulative distribution function of the standard normal distribution. The second variant 
allows to emulate the effect of a non-constant bandwidth for the radial basis function. 

4.2.2 Feasibility Restoration 

Feasible decisions are generated by completing the Gaussian process model by the generic 
projection method fill I) . The program fill I) has no parameter to tune. 

We also study the behavior of a feasibility restoration heuristic well adapted to the test 
problem. The heuristic depends on some priority order over the coordinates j of the vectors 
Qt,j i n ffT4"l) . It consists in creating inventory variables s, and initializing them to the values 
qt-i,i for all i, and then trying to reach the quantities Af ^ of the Gaussian model by consuming 
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the products (t — in the needed proportions. Namely, 



For all i, set Sj = qt-i,i- 

For all j considered sequentially according to a pre-specified order a t , 

define q[~ = mmi{si/ A t>ij : A t;ij > 0}; (17) 

set q t ,j = min{A^,(ft-}; 

and for all i, replace Sj by s« — A t ^q t j. 

The priority orders cr t are viewed as the parameters of the heuristic. We generate the priority 
orders randomly, by sampling in the space of permutations. 

4.3 Scenario- Tree Approximations 

In this section we describe how the scenario trees are built and how the shrinking-horizon 
procedure for out-of-sample validation is implemented. The shrinking-horizon procedure is 
the benchmark against which the learned policies will be compared. 

4.3.1 Method for Constructing the Scenario Trees 

We considered scenario trees with a uniform branching factor b. We used an optimal quan- 
tization approach for choosing the b discrete values for £ t and assigning to them probability 
masses ( jPages and Printemsl . 120031 ) . In a nutshell, this approach works by selecting values 



^ that minimize the quadratic distortion 

^ 2 ({d 4) }ti) = %{min||ef ) -6l| 2 }, 

1<1<0 

and then evaluates probabilities by integrating the density of the distribution of over 
the cells of the Voronoi partition induced by the points £ t . 

A scenario tree on T stages has b T ~ x scenarios (exponential growth). By solving scenario- 
tree approximations on trees with increasing branching factors, we determined that the test 
problem could be solved to a reasonable accuracy with a branching factor b = 10 (Figure |3]). 
The solving time grows exponentially with b. 

4.3.2 Value of the Multistage Model 

The optimal value of the multistage model is about -375 (corresponding to a net profit), cf. 
Figure [3j This value should be compared with the optimal value of the two-stage model ffT5|) . 
Our estimate for it is -260. The value of the multistage model over the two-stage model is 
thus in theory an expected profi t increase of about 45 perce nt. 



We recall that according to iBirge and Louveauxl (Il997l ). there does not seem to be a 



structural property in multistage models that would guarantee a large value over their deter- 
ministic or two-stage counterpart; the value is very dependent on the numerical parameters. 
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Branching Time 

factor b (in seconds) 

3 5^0 

4 8.0 

5 17.2 

6 40.7 

7 79.8 

8 177.5 

9 353.5 
10 670.6 

Figure 3: Scenario-tree approximation: (Left) Convergence of the optimal value to the true 
optimal value, and (Right) solution time. 

4.4 Computational Results for the Benchmark Policies 

Using (1121) on a fixed test set of 10000 scenarios, we have estimated the value of solving a new 
scenario-tree approximation of the problem (fl4"l) at each stage over the remaining horizon, 
given the current information state. To investigate the effect of the size of the scenario tree on 
the policy, we have tested the procedure on 3 choices for the branching factor: 6 = 3,5, 7. The 
first-stage decision is constant and is computed only once (with the chosen value of b). Then 
the simulation of the decision process for optimizing the recourse decisions X2,X3,x^ online 
is run on 10000 scenarios, using a scenario tree with branching b over the remaining horizon. 
Actually for x 4 we use the closed-form solution x^i = min{g 3i j, r/i}. The result of these 
simulations is presented in Figure H] for the performance of the policies (curve "Benchmark") 
and in Tabled] for the computation time on a single processor (row "Benchmark"). 

4.5 Computational Results for Learned Policies 

We have reported on Figure H]the results of simulations on 10000 scenarios for three variants: 

1. GP-1: covariance function f fl6|) with g(z) = z, feasibility restoration procedure (fTTj) . 

2. GP-2: covariance function f )16p with g(z) = $(z), feasibility restoration procedure fllip . 

3. GP-3: covariance function (1161) with g(z) = §(z), feasibility restoration procedure (|T7|) . 

Each variant was tested on the three datasets collecting the state-decision pairs from a 
scenario tree with branching factor 6 = 3,5, 7. The determination of the best algorithm 
parameters for each variant was made by direct search, treating each set of parameters as a 
possible model for the decision policy. 

On the test problem we have considered, it seems that the simple program (TTTT) introduces 
a large computational overhead. The simulation times of the models GP-1 and GP-2 are 
only 1.5 to 6 times faster than the benchmark method. One possible explanation is that 
the scenario-tree approximations built on the remaining horizon have not many scenarios 
and thus are not really much more difficult to solve than the myopic program (TTTT) . When 
we replace (ITT|) by the heuristic (IT7I) in GP-3, we obtain a very important speed-up of the 



-375 - 
-385 
-395 - 




3456789 10 
Branching factor b 
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-345 
-355 
-365 
-375 



-K . 



. -x 



GP-1 
GP-2 

GP-3 

Benchmark 



3 5 7 

Figure 4: Simulation of policies on 10000 scenarios: Averaged value 



Table 1: Simulation of policies on 10000 scenarios: Computation time 



Policy 


Cpu time (in seconds) 




6=3 6=5 6=7 


Benchmark 


17000 31000 65000 


GP-1 


11000 11000 12000 


GP-2 


12000 12000 12000 


GP-3 


10 10 10 



simulations (Tabled]) for a relatively small loss of performance with respect to the benchmark 
(Figure H]) and the theoretical best value of the multistage program (Figure [3] with b = 10). 

From these experiments, we conclude that the most attractive forms of learned policies 
(decision rules) for doing simulations on a large number of scenarios are those that elimi- 
nate completely the calls to a solver, even for solving very simple programs. Doing so can 
be achieved by setting up a feasibility restoration heuristic that need not be very clever, 
because the implemented decisions depend mostly on the predictions obtained from the 
nonparametric Gaussian process model. 



5 Application: Optimal Selection of Scenario Trees 

This section presents our solution approach to multistage problems over long horizons, based 
on the selection of random scenario trees by out-of-sample validation. Section 15.11 makes 
the stylized example of Section [2] more concrete. Section 15.21 presents an algorithm for 
generating random branching structures, that turned out to be well adapted to the problem 
of Section 15. II Sec t ion 15.31 presents the computational study. This material is also presented 



m 



Defourny et al.l (120121 ) . Section |5~41 considers the number of random trees to sample. 



16 



5.1 Studied Problem 



Consi der the following problem, studied by lHilli and Pennanenl (120081 ) and lKiichler and Vigerske 



( 120101 ). where p is a risk-aversion parameter and r\ a budget parameter: 
min^,...^ p~ 1 logE{exp{-p^^ =1 ^ x t (£ [t] )}} 

subject to ac t (f[ t ]) G «*i(o;i(^[i]), . . . , a; t _i(^ t _i])), t = 1, . . . ,T, (18) 
Xt(xi(£[i]),...,x t -i(€[t-i])) = { x t eR:J2T=i x r(^[T])+x t <rj, < x t < 1 } . 

The random variables £i, • • • , £t are generated from 

& = s t - k, St = s t _i exp{(xei - cr 2 /2} with s = k , (19) 

where a 2 = 0.07, k — 1, and ei,...,e 4 are i.i.d. from the standard normal distribution. 
(Since £i is truly stochastic, the problem is over T recourse stages. We could introduce a 
trivial constant first-stage decision xq = associated to £o = 0, so strictly speaking the 
multistage model is over T + 1 stages.) 

When p tends to 0, the program ([TBI) becomes linear, and for this case an optimal policy 
is the simple bang-bang policy: Xf = 1 if £ t > and t > T — rj; x t = otherwise. 

Our goal is to solve (1181) on T = 52 for various values of p and r]. On long horizons, it is 
out of the question to consider scenario trees with uniform branching factors (with b = 2 we 
have already 2 52 = 4.5 • 10 15 scenarios). 

Interestingly, the decisions optimal for scenario-tree approximations of ([TBI turn out to 
be very dependent on the branching structure of the tree. When a branching is missing in 
one scenario of the tree, a deterministic vision of the future is induced for that scenario from 
the stage of the missing branching to the stage where a branching is finally encountered. This 
will not hurt if the value of the multistage model on this part of the scenario and onwards 
is negligible, but we cannot know that in advance (that is, prior to having computed an 
optimal policy that solves (ITS]) or at least having determined its structure). 

There does not seem to be much advantage in devoting computational resources to an 
optimization of the branching structure of the tree, because at the end of the day we would 
still be unable to estimate how realistic the optimal value of the approximation is with 
respect to the true optimal value or with respect to a binary tree of 2 52 scenarios. 

Motivated by these considerations, we propose to generate branching structures purely 
randomly. The approach makes sense only if we can estimate the value of the approximation 
for the true problem without an optimistic bias. This is in turn possible by simulating 
a decision policy learned from the solution to the scenario-tree approximation. With a 
learning and simulation procedure that is fast enough, we can score several scenario trees 
and thus in essence implement a Monte Carlo search algorithm in the high- dimensional space 
of branching structures. 



5.2 Random Generation of Branching Structures 

Based on various numerical experiments for solving ([TBI) , the algorithm described in Figure [5] 
has been found to work well for generating, with a sufficient probability, branching structures 
leading ultimately to good decision policies. 
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Given N (desired approximate total number of scenarios): 

1. Create a root node (depth 0). Set t — 0. 

2. Set v t to the number of nodes at depth t. Set r t = (l/u t )(N — 1)/T. 

3. For each node j of depth t: 

Draw Zj uniformly in the interval [0,1]. 

If Zj < r t , append 2 children nodes to node j (binary branching). 
If Zj > r t , append 1 child node to node j (no branching). 

4. If t < T — 1, increment t and go to Step 2. 
Otherwise, return the branching structure. 



Figure 5: Random generation of scenario tree branching structures 



The algorithm produces branching structures leading to trees having approximately 
iV scenarios in the following sense. Assume that the number Vt-i of existing nodes at 
depth T — 1 is large. From each node, create one or two successor nodes randomly (refer 
to Step 3 in Figure [5]). By the independence of the random variables Zj that determine the 
creation of 1 or 2 successor nodes, and by the weak law of large numbers, the created random 
number of nodes at depth T is approximately equal to 

u T = z/ T _ x (2 ■ rt-i + 1 ■ (1 - rvi)) = i/ T -i(l + r t _i) 
= v T ^(l + [l/vr-x){N - 1)/T) = z/ T _x + (N- 1)/T. 

Iterating this recursion yields = + T(N — 1)/T = N. To establish the result, we have 
neglected the fact that when v t _\ is small, the random value of v t conditionally to v t -i should 
not be approximated by the conditional mean of v t) as done in the recursive formula. The 
error affects mostly the first depth levels of the tree under development. We have found, by 
generating random trees and estimating the expectation and variance of the number of leafs, 
that the error had a small effect in practice. 



5.3 Computational Results 

We have considered three sets of 25 scenario trees generated randomly over the horizon 
T = 52: the first set with N = T, the second set with iV = 5T, and the third set with 
iV = 25T. The random structures are generated by the algorithm of Section 15.21 The 
scenarios use values of £ t generated randomly according to f|T9|) . 

The inference of decision policies from the solution of the scenario-tree approximations 
uses a more compact representation of the information state t han More de tails on the 



inference and feasibility restoration procedures can be found in iDefournyl (120 10f ) . 

The computational results are summarized in Table [2] for the accuracy (the best values are 
indicated in bold), and in Table [3] for the overall computational complexity of the approach, 
that involves generating the 25 random trees, solving them, and simulating 5 candidate 
policies per tree on 10000 new scenarios. The reported times are relative to a Matlab 
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Table 2: Value of the best policies found for instances of ( TPS!) with T = 52. 
Problem Simulation on 10000 new scenarios: Average value 



p 7] Benchmark Best learned policies, for 3 tree sizes 









N = 52 


N = 260 


N = 1300 


2 


-0.40 




-0.34 


-0.32 


-0.39 


6 


-1.19 


a 
+^ 


-1.07 


-1.03 


-1.18 


20 


-3.64 


a 

o 


-3.59 


-3.50 


-3.50 


0.25 2 


-0.34 




-0.32 


-0.31 


-0.33 


6 


-0.75 




-0.78 


-0.78 


-0.80 


20 


-1.46 




-1.89 


-1.93 


-1.91 


1 2 


-0.22 




-0.25 


-0.22 


-0.24 


6 


-0.37 




-0.53 


-0.53 


-0.54 


20 


-0.57 




-0.96 


-0.98 


-0.96 



Table 3: Computation times 



Problem 


Total cpu time (in 


seconds) 


P V 


N = 52 


N = 260 


N = 1300 


2 


415 


551 


1282 


6 


435 


590 


1690 


20 


465 


666 


1783 


0.25 2 


460 


780 


2955 


6 


504 


1002 


4702 


20 


524 


1084 


5144 


1 2 


485 


986 


4425 


6 


524 


1095 


5312 


20 


543 


1234 


6613 



implementation, run on a single processor, but the nature of our randomized approach 
makes it very easy to parallelize. Our benchmark is the simulation of the bang-bang policy 
optimal for the risk-neutral case p = 0. 

It is somewhat surprising to see that multiplying the number of scenarios by 25 does not 
translate to significantly better results, as shown by comparing the column iV = 52 to the 
column N = 1300 in Table [2j Note, however, that the results with N = 52 are obtained for 
a particular tree of the set of 25. Most of the time, the results on trees with N = 52 are 
poor. Also, having 52 scenarios or 1300 in the tree is equally terribly small compared to the 
exponential number required to solve with certainty the program on T = 52. 
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5.4 On the Required Number of Random Scenario Trees 



Finally, we note that if there exists a randomized algorithm able to generate with probability 
g > a tree from which a "good" decision policy can be learned (we discuss the sense of 
"good" below), then the number M of trees that have to be generated independently for 
ensuring with probability 5 that at least one of them is "good" is equal to 



M{5) 



log(l-*) 
log(l - g) 



(20) 



we 



For instance, if the randomized algorithm generates a good tree with probability 0.01 
need a set of 300 random trees to obtain a good one with probability 0.95. 

The sense of "good" can be made precise in several ways: by defining an aspiration 
level with respect to a lower bound on the tru e valu e of the multistage program, obtained 
for instance with the techniques of iMak et al.l (Il999l ); by defining an aspiration level with 
respect to a benchmark solution that the decision maker tries to improve; or by defining 
aspiration levels with respect to risk measures besides the expectation. 

Indeed, it is possible to compare policies on the basis of the empirical distribution of their 
cumulated cost on a large test sample of independent scenarios. 



6 Conclusions 

This paper has presented an approach for inferring decision policies (decision rules) from the 
solution of scenario-tree approximations to multistage stochastic programs. Precise choices 
for implementing the approach have been presented in a Bayesian framework, leading to a 
nonparametric approach based on Gaussian processes. The sensitivity of the approach has 
been investigated on a particular problem through computational experiments. 

The inference of decision policies could be a useful tool to calibrate scenario tree gener- 
ation algorithms. This line of research has been followed by developing a solution strategy 
that works by generating scenario trees randomly and then ranking them using the best 
policy that can be inferred from their solution. Further work could be useful for identifying 
randomized algorithms likely to generate good scenario trees. If these algorithms exist, a 
solution strategy based on them could fully leverage the computing power of current super- 
computer architectures. 
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A Numerical Parameters 



The value of the numerical parameters in the test problem (JT4"j) are given here. 



ci 

C2 
C3 
C 4 
h 
h 
h 

w 
h 



A, 



[ 0.25 1.363 0.8093 0.7284 0.25 0.535 
[ 2.5 2.5 2.5 2.5 13.22 2.5 3.904 2.5 
[ 3.255 2.5 2.5 8.418 2.5 ] T x el 5 
- [ 21.87 98.16 31.99 10 10 ] T 
[ 13.9 9.708 2.14 4.12 ] T 



0.25 

T 



0.25 
x G I 



0.25 0.4484 

d8 ^/ tu>12x8 



0.25 0.25 



p8x5 



12.86 


9.9U1 


6.435 


7.446 


1 

J 








18.21 


7.889 


3.2 2.679 ] 


r 








10.14 


4.387 


9.601 


4.399 


] 








17.21 


4.983 


7.266 


9.334 










0.4572 





4.048 











0.8243 


11.37 








0.7674 


0.5473 


0.3776 











0.4794 





0.4861 


1.223 





1.475 




















0.5114 


0.3139 











12.29 


1.378 





0.3748 


0.4554 








0.7878 





0.293 


1.721 














1.504 


0.4696 


0.248 





0.1852 





0.3486 








1.204 





0.7598 


0.452 

















0.2515 


0.3753 


0.6249 





1.248 





1.545 

















0.2732 














0.6597 





2.525 














1.595 








1.51 


1.041 


0.9847 





1.223 


0.6367 





o - 

















1.111 

















0.4579 

















0.1693 


0.6589 














0.5085 


2.643 

















0.4017 























0.7852 


85.48 























0.806 


0.5825 . 
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